!>\brief
!! This is a simple test case on the square \f$[0\,,2\pi]^d\f$, using
!! periodic boundary conditions.
!!
!! \n
!!
!! We consider the problem
!! \f{displaymath}{
!!  \begin{array}{rcll}
!!   \Delta u & = & f & {\rm
!!   in}\,\,\,\Omega
!!  \end{array}
!! \f}
!! with \f$\Omega=[0\,,2\pi]^d\f$ and periodic boundary conditions.
!! The exact solution is then
!! \f{displaymath}{
!!  u(x) = \prod_{i=1}^d \cos\left(ix_i+i\right),
!! \f}
!! with gradient
!! \f{displaymath}{
!!  \partial_{x_i}u(x) = -i\sin\left(ix_i+i\right)
!!  \prod_{\begin{array}{c}j=1\\j\neq i\end{array}}^d
!!  \cos\left(jx_j+j\right),
!! \f}
!! total flux
!! \f{displaymath}{
!!  q_i(x) = \partial_{x_i}u(x)
!! \f}
!! and divergence of the total flux, equal to the forcing term,
!! \f{displaymath}{
!!  \nabla\cdot q(x) = f(x) = -\left(\sum_{i=1}^di^2\right)u(x).
!! \f}
!<----------------------------------------------------------------------
module mod_periodic_test

!-----------------------------------------------------------------------

 use mod_messages, only: &
   mod_messages_initialized, &
   error,   &
   warning, &
   info

 use mod_kinds, only: &
   mod_kinds_initialized, &
   wp

!-----------------------------------------------------------------------
 
 implicit none

!-----------------------------------------------------------------------

! Module interface

 public :: &
   mod_periodic_test_constructor, &
   mod_periodic_test_destructor,  &
   mod_periodic_test_initialized, &
   test_name, &
   test_description,&
   coeff_diff,&
   coeff_adv, &
   coeff_re,  &
   coeff_f,   &
   coeff_lam, &
   coeff_q,   &
   coeff_divq

 private

!-----------------------------------------------------------------------

! Module types and parameters

 ! public members
 character(len=*), parameter :: &
   test_name = "periodic"

! Module variables

 ! public members
 character(len=100), protected :: &
   test_description(1)
 logical, protected :: &
   mod_periodic_test_initialized = .false.
 character(len=*), parameter :: &
   this_mod_name = 'mod_periodic_test'

!-----------------------------------------------------------------------

contains

!-----------------------------------------------------------------------

 subroutine mod_periodic_test_constructor(d)
  integer, intent(in) :: d

  character(len=*), parameter :: &
    this_sub_name = 'constructor'

   !Consistency checks ---------------------------
   if( (mod_messages_initialized.eqv..false.) .or. &
          (mod_kinds_initialized.eqv..false.) ) then
     call error(this_sub_name,this_mod_name, &
                'Not all the required modules are initialized.')
   endif
   if(mod_periodic_test_initialized.eqv..true.) then
     call warning(this_sub_name,this_mod_name, &
                  'Module is already initialized.')
   endif
   !----------------------------------------------

   write(test_description(1),'(a,i1)') &
     'Periodic test case, dimension d = ',d

   mod_periodic_test_initialized = .true.
 end subroutine mod_periodic_test_constructor

!-----------------------------------------------------------------------
 
 subroutine mod_periodic_test_destructor()
  character(len=*), parameter :: &
    this_sub_name = 'destructor'
   
   !Consistency checks ---------------------------
   if(mod_periodic_test_initialized.eqv..false.) then
     call error(this_sub_name,this_mod_name, &
                'This module is not initialized.')
   endif
   !----------------------------------------------

   mod_periodic_test_initialized = .false.
 end subroutine mod_periodic_test_destructor

!-----------------------------------------------------------------------
 
 pure function coeff_diff(x) result(mu)
  real(wp), intent(in) :: x(:,:)
  real(wp) :: mu(size(x,1),size(x,1),size(x,2))

  integer :: i
 
   mu = 0.0_wp
   do i=1,size(x,1)
     mu(i,i,:) = 1.0_wp
   enddo
 end function coeff_diff
 
!-----------------------------------------------------------------------
 
 pure function coeff_adv(x) result(a)
  real(wp), intent(in) :: x(:,:)
  real(wp) :: a(size(x,1),size(x,2))

   a = 0.0_wp
 end function coeff_adv
 
!-----------------------------------------------------------------------

 pure function coeff_re(x) result(sigma)
  real(wp), intent(in) :: x(:,:)
  real(wp) :: sigma(size(x,2))
 
   sigma = 0.0_wp
 end function coeff_re
 
!-----------------------------------------------------------------------
 
 pure function coeff_f(x) result(f)
  real(wp), intent(in) :: x(:,:)
  real(wp) :: f(size(x,2))

  integer :: i

   f = real(sum( (/(i , i=1,size(x,1))/)**2 ),wp)*coeff_lam(x)

 end function coeff_f
 
!-----------------------------------------------------------------------

 pure function coeff_lam(x) result(l)
  real(wp), intent(in) :: x(:,:)
  real(wp) :: l(size(x,2))

  integer :: i
  real(wp) :: ri

   l = 1.0_wp
   do i=1,size(x,1)
     ri = real(i,wp)
     l = l * cos(ri*x(i,:)+ri)
   enddo
 end function coeff_lam
 
!-----------------------------------------------------------------------

 pure function coeff_q(x) result(q)
  real(wp), intent(in) :: x(:,:)
  real(wp) :: q(size(x,1),size(x,2))

  integer :: i, j
  real(wp) :: ri, cx(size(x,1),size(x,2)), sx(size(x,1),size(x,2))

   do i=1,size(x,1)
     ri = real(i,wp)
     cx(i,:) = cos(ri*x(i,:)+ri)
     sx(i,:) = sin(ri*x(i,:)+ri)
   enddo
 
   do i=1,size(x,1)
     q(i,:) = real(i,wp)
     do j=1,size(x,1)
       if(j.eq.i) then
         q(i,:) = q(i,:) * sx(j,:)
       else
         q(i,:) = q(i,:) * cx(j,:)
       endif
     enddo
   enddo
 
 end function coeff_q
 
!-----------------------------------------------------------------------

 pure function coeff_divq(x) result(dq)
  real(wp), intent(in) :: x(:,:)
  real(wp) :: dq(size(x,2))

   dq = coeff_f(x) 

 end function coeff_divq
 
!-----------------------------------------------------------------------
 
end module mod_periodic_test

